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Abstract 



We use a random matrix model to study chiral symmetry breaking in QCD at 
finite chemical potential /i. We solve the model and compute the eigenvalue 
density of the Dirac matrix on a complex plane. A naive "replica trick" 
fails for /i 7^ 0: we find that quenched QCD is not a simple n ^ limit of 
QCD with n quarks. It is the limit of a theory with 2n quarks: n quarks 
with original action and n quarks with conjugate action. The results agree 
with earlier studies of lattice QCD at 7^ and provide a simple analytical 
explanation of a long-standing puzzle. 
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I. INTRODUCTION 



The spontaneous breaking of chiral symmetry is one of the most important dynamical 
properties of QCD which shapes the hadronic spectrum. A great deal of understanding 
of this nonperturbative phenomenon at zero and finite temperature has been achieved by 
various methods In particular, we expect that the chiral symmetry is restored above a 
certain critical temperature. The study of this new chirally symmetric phase of hot QCD 
is one of the primary objectives of heavy ion colliders. In contrast, the behavior of QCD 
at large baryon density (conditions which can arise in the heavy ion colliders or in neutron 
stars) is not well understood. The main puzzle has for a long time been a contradiction 
between a straightforward physical expectation and numerical results from quenched lattice 
QCD Simulations with dynamical quarks, on the other hand, are very inefficient at 

finite fi — the fermion determinant is complex. 

The puzzle concerns the dependence of the order parameter (the chiral condensate {ipip)) 
on the baryon chemical potential. A non-analytical change in the value of {ipip) should occur 
when yU > ~ rriB/S, where tub is the mass of the lightest baryon. At this point the 
production of baryons becomes energetically favorable. For smaller fi the value of (ipip) is 
nonzero. In contrast, lattice simulations of quenched QCD indicate that yUc = (at zero 
bare quark mass), i.e., the chiral condensate vanishes if yU 7^ |@J^. A number of possible 
explanations has been suggested However, the answer to this puzzle remains unclear. 

This work was motivated by a desire to shed some light on this question using the 
random matrix approach which received considerable interest recently p|-p!2|. It is based 
on the idea that, for the purpose of studying chiral symmetry breaking, fluctuations of the 
Dirac operator in the background of the gauge flelds can be approximated by purely random 
fluctuations of its matrix elements in a suitable basis. For example, in the instanton liquid 
model this basis can be formed from the Dirac zero modes for individual (anti)instantons, 
which due to overlaps form a band of small eigenvalues responsible for the chiral symmetry 
breaking |jl3|. A similar random matrix approach is fruitful in the studies of spectra of 



systems with a high level of disorder, such as spectra of heavy nuclei [|14|. Introduction of 
chemical potential into such a model of chiral symmetry breaking is straightforward. The 
resulting Dirac matrix (times i) is non-hermitian. Thus the eigenvalues lie in the complex 
plane rather than on a line. Such random matrix models have not received much attention 
previously and this study is a step in an unexplored direction. 

In this Letter we show how to solve such a model in the thermodynamic limit and discuss 
the implications. 
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II. THE RESOLVENT 



In order to study chiral symmetry breaking we shall calculate the resolvent of the Dirac 
operator D: 

G=(tTiz-D)-'), (1) 



as a function of the bare quark mass z which we take to be a complex variable z = x + iy. 
The average is over fluctuations of the random matrix elements of -D. It should be obvious 
that G is the same as {'ip'ip) . The resolvent can be expressed through the average eigenvalue 
density p: 

G{x,y) = J dx'dy'p{x',y')^^. (2) 

A vector G = (ReG, — ImG) is the electric field created by the charge distribution p. This 
makes the inversion of (|^) obvious: 

1 - - Id 

p=^G = -—G, (3) 
ivr TT oz* 

where d/dz* = (d/dx + id/dy)/2. 

Analytical properties of G are very closely related to the chiral symmetry breaking. From 
(§) we see that p vanishes if the function G is holomorphic. A discontinuity of G along a 
cut going through 2; = is the signature of the spontaneous chiral symmetry breaking: 
('?/''?/') (+0) 7^ {tptp){—0). This observation together with leads to the Banks-Casher 
relation [0: (ipip) = vrp(O), where p(0) is the density per length on the cut at 2; = 0. 
However, (|^) is more general and can be applied to a case when the non-analyticity is not 
in the form of a cut but occupies a 2-dimensional patch, which is the case in our model. 

III. THE MATRIX MODEL AND NAIVE REPLICA TRICK 

The matrix D has the form: 




where we added the chemical potential term p'-^Q to the Dirac matrix [|5[. The N x N 
matrix elements of X are independently distributed complex Gaussian random variables: 
P(X) = const X exp{— Tr XX"^}. The unit of mass in the model is set by n^/ {^ip)o ~ 
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200 MeV, where 77,4 ~ 1 fm is the number of small eigenvalues in a unit volume (instanton 
density [|l^]) and {ipip)^ ~ (200 MeV)^ is the chiral condensate at T = 0, /i = 0. 



In order to find the resolvent (|l|) we introduce n quark fields (replicas) and calculate: 

K = -^ln(det"(2-D)). (5) 
This quantity continued to n — > (quenched limit) becomes: 

V = - {\n<lei{z - D)) , (6) 

from which we find G: 

""-If- « 

The trace in (|l]) is normalized as tr 1 = Tr 1/ (2A^) . Following the electrostatic analogy of 
the previous section one can view YieV as the scalar potential for G. 
Using Hubbard-Stratonovitch transformation we obtain: 

exp{-nK} = / Vadei^ Iz + a ^l \ ^^t|^ (g) 

\ (1 z + a) j 

where a is an auxiliary complex nxn matrix field. For large the calculation of the integral 
amounts to finding its saddle point. If we assume that the replica symmetry is not broken 
(i.e., a is proportional to a unit matrix) we arrive at the saddle point equation: 

{z + a) = a[{z + af - n\ (9) 

The complex value of a in (|^) is the analytical continuation of the real part of the diagonal 
matrix elements of a in (8). The imaginary part is zero in the saddle point. The solution of 
this cubic equation is straightforward. It is the same as in a similar model |T(]| with uj ifi. 



Finally, it is easy to find using (|^J|) that G = a where a is the saddle point given by (§). 

The Vn does not depend on n and the limit — > seems obvious. However, in the next 
section we shall compare this expectation to numerical data and see that the limit n — is 
in fact very different! Now let us summarize the properties of this model for n > 0. 

We see that G{z) is a holomorphic function. It has 3 Riemann sheets and we select the 
one where G ^ 1/z for z ^ 00 - which follows from / dxdy p = 1 and the Gauss theorem. 
The only singularities on the physical sheet are the pole at z = 00 and 2 (for < 1/8) or 
4 (for > 1/8) branch points connected by cuts. The brunch points are where 2 of the 3 
solutions of (^ coincide. The trajectory of a cut is determined by a condition that a is the 
deepest minimum (out of 3) of: Re [a^ — \n{{z + a)^ — /i^)]. 
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At /i = the cut along imaginary axis connects two singularities at 2; = ±21 For nonzero 
/i the singularities start moving towards each other along the imaginary axis. At = 1/8 
each of the branch points bifurcates in two ones which move off the y axis into the complex 
plane. The cut goes through the origin (along the y axis) until /i^ = 0.278.... At this point 
it splits into two cuts connecting complex conjugate points. This means that //^ = 0.278... 
in such a model. 



IV. NUMERICAL RESULTS AND THE SOLUTION OF THE MODEL 



For n = one can easily determine the density of eigenvalues numerically by calculating 
the eigenvalues of the random matrix D and plotting them on a complex plane. The density 
of points on such a scatter plot is proportional to p. The results for different values of p? are 
shown in Fig. 1. They contradict naive expectations from the previous section. At /i = 
all eigenvalues are distributed between points z = ±2i on the y axis. However, already at 
very small nonzero /x^ <^ 1/8 the eigenvalue density is nonzero in a "blob" of finite width 
in X direction which grows with fi. The same behavior is seen in quenched lattice QCD p| 
and gives rise to the paradox described in the Introduction: there is no discontinuity in the 
value of (ipip) at any fi > 0. The matrix model has an advantage: it is amenable to exact 
treatment which clarifies the nature of the problem. 

The failure of the naive replica approach can be understood if we look at the expres- 
sion (|): it does not contain z*\ On the other hand, eq. (Q) tells us that p 7^ if G depends 
on z*, i.e., if it is not holomorphic. In fact, the correct replica trick for a non-hermitian 
matrix should start from the quantity: 

Vn,n = -- In (det"(z - D){z* - D^)) , (10) 

which is now real due to introduction of the quarks with conjugate Dirac matrix. Naively, 
in the limit n —>■ the conjugate quarks decouple but, as we shall see, this is not always the 
case! In mathematics an analogous construction is called a V-transform [|l^ and allows one 
to study spectra of non-hermitian matrices. In the present context this formal construction 
has a clear and simple physical meaning. 

We can calculate ( [T0| ) using the same method as for (^. Now, however, we have to 
introduce 4 auxiliary complex n x n fields, and we arrive at: 
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exp{-nVn,n) = J VaVhVcVd det^ 



z + a fi 
fi z + ic 
id^ z* + 



id 




2C 







xexp{-iV(|a|2 + |6|2 + |c|2 + |d|2)}. 



(11) 



The set of solutions of the saddle point equation is richer in this case. There is a 
solution with c = d = 0. In this case the conjugate quarks do decouple and we obtain the 
same holomorphic function G as before. However, there is another solution in which the 
condensates c and d are not zero! Then the function G is not holomorphic and therefore 
p 7^ 0. This saddle point dominates the integral at small z for < /i < 1. 

The condensates c and d are bilinears of the type (V'x); mixing original ip and conjugate 
X quarks. These condensates do not break the original chiral symmetry but a spurious 
(replica type) symmetry involving both original and conjugate quarks. Similar condensates 
carrying baryon number were discussed in the SU{2) model of QCD with quarks [|l^. In the 
quenched theory, as in the original chiral symmetry is always restored at /i > 0. The 
spurious symmetry is spontaneously broken for /i < 1 and is restored for p > 1. 

The boundary of the p 7^ region is given by: 



y 



- p') - (1 + 4p'^ - 8fi^)x' - Afi^x^' 



(12) 



It is plotted on Fig. 1 for comparison with numerical data. The baryonic condensates c and 
d inside of the "blob" are given by: 



MP 



p^ 



p^ — 



p 



X 



4(p2 



X 



2^2 



r 
4 



(13) 



On the boundary (|T^ they vanish and the two solutions (holomorphic and non-holomorphic) 
match. In the outer region: c = d = and G = a is the solution of the cubic equation (^. 
Inside of the "blob" the resolvent is given by: 



G 



1 

a = —- 



X 



2p2 



X 



2 ^ 2 ' 



and the density of the eigenvalues (H) is: 



+ p^ 



X 



2^2 



(14) 



(15) 



^ 47r V(p2 

To appreciate non-triviality of this result one should notice that expression ([^) which 
defines the resolvent appears to depend only on z\ The limit — > must be taken with 
great care, as is well-known in the replica approach |]19 . 
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V. CONCLUSIONS 



The fermion determinant in QCD is complex at nonzero chemical potential. Lattice 
simulations of such a theory are extremely inefficient. Therefore all reliable data from lattice 
QCD so far have been obtained for a quenched theory. We learn from the random matrix 
model that the quenched theory at finite fi behaves qualitatively different from the QCD 
with dynamical quarks. Rather, the quenched approximation describes a theory where each 
of the quarks has a conjugate partner, so that the fermion determinant is non-negative. We 
see that for such a theory the result /Xc = is natural. Similar arguments have been given 
by several authors in different settings and using less realistic models [Q. Here it can be 
demonstrated in a very clean and explicit way. 

Simulations with dynamical quarks at strong coupling are possible in SU{2) and SU{4) 
QCD [|18|] and also agree with our results. The fic is finite in the SU{A) theory. On the other 
hand, in the SU (2) theory, where the quarks are self-conjugate, /Xc = due to the baryonic 
condensates. 

The matrix model describes many features of the chiral symmetry breaking in QCD very 
well [^0. One of the apparent limitations, however, is that it is static — there are no 
kinetic terms and we cannot study spectrum of masses. In the quenched QCD fic appears to 
coincide with half of the mass of the so-called baryonic pion 0] — a bound state of a quark 
and a conjugate antiquark. It is degenerate with the tt- meson but carries a nonzero baryon 



number. From the exact solution (|12]) we find /ic ~ ym/2, for small quark mass m <^ 1. If 
we had /tt in our model we could relate fic to the mass of the pion. The model also does not 
account for the confinement of quarks. It remains to be seen if the confinement plays a role 
in the case under consideration. 
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FIG. 1. Scatter plots on a complex plane (x, y) of the eigenvalues of an ensemble of 20 random 
100 X 100 matrices D at 4 values of /u^: 0.06 (a), 0.10 (b), 0.40 (c) and 1.20 (d). The solid curves 
follow eq. (|12|). 
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